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We have studied the collective behavior of a population of integrate-and-fire oscillators. We show 
that diversity, introduced in terms of a random distribution of natural periods, is the mechanism that 
permits to observe self-organized criticality (SOC) in the long time regime. As diversity increases 
the system undergoes several transitions from a supercritical regime to a subcritical one, crossing 
the SOC region. Although there are resemblances with percolation, we give proofs that criticality 
takes place for a wide range of values of the control parameter instead of a single value. 

PACS numbers: 05.90.+m,87.10.+e,64.60.Ht 
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In spite of the great interest received along the last 
decade by many systems exhibiting self-organized criti- 
cality (SOC) (T), it is still an open question to find neces- 
sary or sufficient conditions to observe this phenomenon, 
since there is no framework to predict, a priori, whether 
an arbitrary extended system will be, by its own dynam- 
ics and without any parameter tuning, critical in the long 
time regime. 

Nevertheless, there is a set of common trends which 
characterize systems displaying SOC 0. One of them 
concerns the dynamics that drives the elements of the 
system to a certain threshold. When some unit reaches 
the threshold, interaction between elements takes place, 
triggering a chain process or avalanche that ends when 
all the elements are below the threshold again. Then, 
a power-law distribution of avalanche sizes is the hall- 
mark of SOC. A key point, crucial to observe SOC, is the 
separation between the slow time scale associated with 
the process that leads the units to the threshold (driv- 
ing), and the fast time scale associated with the interac- 
tion (avalanches). Conservation was also believed essen- 
tial to obtain SOC. Certainly, for the sandpile model |l]] 
and other randomly driven models it is an indispensable 
requirement. A nonconservative dynamics introduces a 
characteristic length independent of system size || . How- 
ever, several continuously driven models proposed later 
changed the widespread belief Q. In these models SOC 
is not necessarily destroyed in a nonconservative regime 
and the distribution of avalanche sizes follows a power- 
law decay in a wide region of parameter space, with ex- 
ponents depending on the level of dissipation. 

Another point not studied so profoundly concerns the 
individual features of each element in the system. Up to 
now, it has been common to assume that all the units are 
identical. However, if SOC should have any relevance in 
physics or biology it should be robust in spite of the in- 
herent differences between the members of a population. 
In other words, diversity should not destroy the critical 
properties of a given self-organized system. The object 
of this paper is to show that diversity not only does not 



break SOC but as a matter of fact, it is the mechanism 
that enables to observe it for certain continuously driven 
models which do not exhibit SOC under normal circum- 
stances. The subject has a clear general interest. There 
are some collective phenomena such as the mutual syn- 
chronization of the members of a biological population 
H which traditionally have been tackled by assuming 
that all the units are identical. However, this assump- 
tion is not a necessary requirement. Several authors [^[0] 
have shown that after a suitable modeling, a group of 
nonidentical oscillators, each endowed with its own nat- 
ural frequency, picked from a random distribution, may 
display a coherent temporal activity if the disorder level 
is below a certain critical value. A less intuitive oppo- 
site behavior has been also reported: disorder (diversity) 
can remove chaos and foster synchronization in a certain 
model of oscillators (§) . Uncorrelated differences between 
the members of the population trigger regular spatiotem- 
poral patterns. In this paper we give evidence of another 
related phenomenon. A group of pulse-coupled oscilla- 
tors evolve in a complex manner if they are identical, 
generating avalanches with many characteristic sizes (re- 
lated with the linear dimension of the system). However, 
diversity will change the collective properties of the long 
time regime and will induce SOC. 

Let us consider a population of integrate-and-fire os- 
cillators. Each oscillator is defined in terms of a state 
variable E which evolves in time as 



dEi 



(i) 



and when Ei reaches a threshold value E t h, the i-th os- 
cillator relaxes, and Ei is redistributed instantaneously 
among its neighbors (labeled by n) according to 



Ei 
E n 
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and so on for every Ek > E t h Vfc. This process, that 
continues until Ek < Eth Vfc, constitutes an avalanche 
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whose size s is given by the number of relaxations . If 
ei n < E t h^i except for at least one element for which 
the inequality is strict, it is guaranteed that avalanches 
of infinite size are impossible. Since the relaxing site is 
reset to zero and a fixed quantity Ei n is transferred to 
the neighbors, the model is intrinsically nonconservative. 
It is assumed jEth < S and 7 is a non-negative con- 
stant whose physical meaning depends on the model one 
is dealing with. For instance, this model mimics a sim- 
plified version of the dynamics of spiking neurons, ideal- 
izing the cell membrane as an RC circuit || . Ei denotes 
the membrane potential of a given neuron, 7 -1 = RC 
the membrane time constant, S (in appropriate units) 
a constant current that does act as a driving, and Ei n 
the synaptic coupling strength between neurons i and n. 
Equations (Q) and may also model the evolution of 
the cardiac pacemaker JhJ, swarms of flashing fireflies, 
and many other biological systems [[IT 12|. 

In order to study local connectivity we have consid- 
ered the case of a two-dimensional square lattice of lin- 
ear size L with nearest-neighbor uniform interactions, 
ei n = e. Recent studies on integrate-and-fire neurons 
jl3| are more devoted to lattice models with periodic 
boundary conditions. However, the assumption of open 
boundary conditions breaks the homogeneous connectiv- 
ity allowing the boundary units to be connected with less 
neighbors than the bulk units. This assumption will be 
present in the rest of the paper and makes the model 
also interesting in other fields. For 7 = it reduces to 
a coupled map lattice proposed by Feder and Feder as a 
stick-slip model of earthquakes Q . 

In addition, we have introduced diversity in terms of 
a random distribution of intrinsic periods. The period of 
each oscillator is given by 
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There are different ways to introduce such type of 
quenched disorder in the model. One possibility is to 
assume a distribution of R and C. Another option is to 
consider a distribution of input currents S. Both situa- 
tions are plausible from a realistic point of view, but we 
have considered the latter. Let us mention that the most 
usual way to introduce diversity in this sort of models is 
by assuming a quenched random distribution of thresh- 
olds |l5|,|l6| . Although the distribution of thresholds also 
implies diversity in the intrinsic periods, it has influence 
in both, the slow and the fast time scale, while our ap- 
proach only affects the slow dynamics. This difference is 
crucial as we will see later. 

When all the oscillators are identical the model de- 
scribed by (Q) and (|2|) does not display SOC for any value 
of 7 > and < e < 0.25. Starting with uniformly dis- 
tributed random initial conditions, Ei G [0, Eth = 1], 
for 7 = and integer ratio Eth/e only large avalanches 



take place because many units reach the threshold si- 
multaneously |i~7[|lg| ]. If the ratio E t h/s is not an inte- 
ger, avalanches of all sizes are observed, but they are not 
power-law distributed JlS} ]. For 7 > (convex driving) 
the model exhibits a complex behavior which, depending 
on the particular values of 7 and e ranges from synchro- 
nization (in the sense that all the avalanches are exactly 
of size L 2 ) to events of all sizes distributed in a compli- 
cated way, as Fig. 1 illustrates. Here we observe that 
the probability density P(s) for an avalanche of size s 
presents a series of peaks at positions that are propor- 
tional to the linear size of the system L. This is a clear 
effect of the open boundaries. Moreover the large peak 
of order L 2 confirms the tendency to synchronization for 
7 > 0, which, however, in this case, the system is not 
able to sustain |Q9|1 . 




FIG. 1. Log- log plot of the stationary distribution of 
avalanche sizes P(s) versus s for the model without diversity. 

The situation changes completely for nonidentical os- 
cillators. For simplicity we have considered a uniform 
distribution of periods. The width A, expressed as the 
length of the symmetric interval (T - A/2,T + A/2) 
centered without loss of generality around T = 1, is a 
measure of disorder or diversity. In Fig. 2 we plot the 
distribution of avalanche sizes for different values of A 
for the same 7 and e as in Fig. 1. We observe several 
stages. First of all, the sequence of peaks displayed in 
Fig. 1 typical of identical oscillators continuously dis- 
appears when diversity increases. Then the distribution 
of avalanches becomes smoother, without intermediate 
peaks, but still maintaining the large one corresponding 
to avalanches of almost the size of the system (L 2 ), as 
displayed for A = 0.15 in Fig. 2 where this trend to- 
wards synchronization can be seen clearly. The behavior 
is supercritical, because there are many events able to 
span the system. More interesting transitions take place 
as disorder increases. For a larger width, the system self- 
organizes in a critical state, without any spatial charac- 
teristic scale, as the power-law distribution of avalanche 
sizes in the curve with A = 0.5 of Fig. 2 indicates. The 
effect of the different periods is to reduce the probability 
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of having large avalanches. Then, for very wide distribu- 
tions of periods one could expect a strong decay of P(s). 
In fact, this is what happens when A > 1.5 for the case of 
Fig. 2, as exemplified by the curve with A = 2: a char- 
acteristic scale independent of the system size appears 
and is responsible of the exponential decay. This means 
that for large diversity avalanches are localized and the 
system is subcritical. These transitions are not sharp, 
and the reported values of A can change with e and 7. 
In particular the loss of criticality and the appearance 
of a finite correlation length has been found difficult to 
characterize. Notice the resemblance between the three 
curves in Fig. 2 and those found in percolation, where 
the critical region is restricted to a single point of the 
control parameter [^o|. However our model shows a fi- 
nite region of criticality instead of an infinitesimal one, 
as we are going to show below. This kind of behavior is, 
for the best of our knowledge, the first case where SOC 
is found between a supercritical region and a subcritical 
one in this class of models. 
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FIG. 2. Same as Fig. 1 for different degrees of diversity. 

Let us pay some attention to the region where the 
power-law decay of P(s) is reported. We have performed 
a finite-size scaling analysis for different widths A. The 
results are shown in Fig. 3. A data collapse for differ- 
ent system sizes L is obtained when plotting L /3 P(s, L) 
against the rescaled variable s/L v . The increment of 
system size does not show any deviation from the scal- 
ing for separate enough values of the control parameter 
A, supporting our statement of a critical region instead 
of a critical point. As a complement we plot in Fig. 
4 the mean size of the avalanches as a function of sys- 
tem size, for different values of diversity. The behavior 
< s >~ L 2v ~P (consequence of the scaling ansatz) even 
for large L confirms the scaling in the critical region. 
In addition, we have released the restriction of identical 
coupling strengths, introducing randomness in space (eij 
quenched random variable) or in space and time {£ij{t) 
annealed), by means of a uniform random distribution 
around the mean value e. We have verified that SOC is 



robust under this perturbation and hence identical cou- 
plings are not a necessary condition to obtain criticality. 
This feature could be relevant in realistic models of spik- 
ing neurons. 

The results shown so far are not characteristic of a 
particular value of the parameters which describe the sys- 
tem. In fact there is a region in the (7,£)-space where 
diversity induces SOC and it corresponds to large val- 
ues of e and small 7. It would be very interesting to 
have knowledge of the complete phase diagram of the 
model. However, taking into account that three param- 
eters are involved: 7, s, and the width A, and different 
system sizes are needed, a complete sweep of the phase 
space would require an enormous effort, that is beyond 
our possibilities. Nevertheless, let us mention that for 
large 7 there exists a range of e values which give com- 
plete synchronization, in the sense previously explained, 
no matter the width of the distribution of periods. 
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FIG. 3. Finite-size scaling analysis of the distribution of 
avalanche sizes for the critical region. For A = 0.3 we obtain 
v — 2.0 and j3 = 3.15 (this curve has been shifted one decade 
upwards for clarity sake) whereas for A = 0.6, v = 1.8 and 
/3 = 2.85 are used. 

Our results have also sense in the context of earth- 
quakes if we imagine the Feder and Feder model as 
a rough version of the Burridge-Knopoff spring-block 
model pjjj . The different intrinsic periods of each unit 
will be caused by different elastic constants in the springs 
connecting the blocks with the driving plate. When 7 > 
a nonlinearity in the elastic response of these springs is 
introduced. With the same goal in mind we have ex- 
amined our disordered model replacing @ by the Olami 
et al. (OFC) rules and we have found that the 

SOC region is robust in spite of a very large disorder, al- 
though eventually it can give rise to localized avalanches. 
These results contrast with the studies performed in Refs. 
fll5| , [i~6| where a distribution of thresholds was considered. 
It was found that while disorder is irrelevant in the con- 
servative regime, it destroys criticality for the dissipative 
case, leading to an exponential distribution of avalanche 
sizes fl5| . A similar change in the collective properties of 
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the disordered system has been used to claim the lack of 
robustness of OFC as a model of earthquakes |l6| . Other 
authors p2j have considered the influence of defects in 
the model. The main result was to observe that SOC is 
robust even for a large number of defects. Notice that 
for the same model randomness included in different pa- 
rameters leads to different behaviors. 
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Linear size of the system L 
FIG. 4. Mean size of the avalanches as a function of the 
system size for different values of diversity using 7 = 0.5, 
e = 0.24. For identical units < s > scales as L 2 . This expo- 
nent decreases continuously with increasing diversity in the 
supercritical region. Notice the accumulation of data points 
for a wide range of A values (starting for A > 0.3) in a nar- 
row interval of mean avalanche sizes, corresponding to the 
critical region. Here the data fit the scaling L 2v ~^ . It is diffi- 
cult to precise up to which values the scaling holds, but when 
A > 1.5 the growth of < s > with L is clearly logarithmic, in 
agreement with a subcritical region. 

Finally, let us remark on the effect that different types 
of noise may have on the collective features of the model. 
The original properties of the Feder and Feder model 
(with 7 = 0) are not robust to noise, e.g., altering the 
relaxation rule by adding a small random number to 
any reset unit, changes the cooperative behavior of the 
system p7| , ^8[ . It does not tend to form a few groups 
of elements with the same phase, but it goes towards a 
SOC state. Note that this type of noise has a completely 
different nature than the quenched source of diversity 
considered in this paper. While the first can be triggered 
by internal fluctuations, the second is an inherent feature 
of each member of the population. Furthermore, while 
the dynamic noise only induces SOC in the linear regime 
(7 = 0), and for very small noise intensities, diversity 
induces SOC in a wide region of the parameter space. 

In summary, we give an example showing that diversity 
is a new mechanism for the emergence of SOC and that 
criticality in noncquilibrium systems is not just a singu- 
larity in parameter space, as it happens in equilibrium. 
Our results have interest for models of integrate-and-firc 
neurons as well as for earthquakes. 
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Note: Just about submitting this paper we became 
aware of Rcf. Q where disorder is introduced in the 
couplings for the OFC model (with 7 = 0), attaining a 
collective behavior in agreement with our results. 
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